function [E e] = xrgr (y, x, mode, Cy)

   ## usage:  [E e] = xrgr (y, x, mode, Cy)
   ##
   ## expanded regression of y on x, optionally returning the error

   if ((nx=columns(x)) < (ny=columns(y)))
      error("xrgr: x dimension too low: x = %d < y = %d", nx, ny) ;
   endif

   global PAR

   E = nan(columns(x), columns(y)) ;
   Jx = any(!isnan(x)) ; x = x(:,Jx) ;

   y = center(y) ; x = center(x) ;

   if (exist("PAR", "var") && isfield(PAR, "wgt") && !isempty(PAR.wgt))
      Wgt = diag(PAR.wgt) ;
   else
      Wgt = eye(columns(y)) ;
   endif 
   y = y * Wgt ;

   if (nargin < 3), mode = "xds"; endif 

   if (mode == "xds")

      w = covest(x) ;
      if any(isnan(w(:)))
	 error("too few x values") ;
      endif 
      w = (w' + w)/2 ;
      GxI = chol(w)^(-1) ;
      if (nargin < 4)
	 w = covest(y', y) ;
      else
	 ##I = all(isnan(y)) ;
	 ##w = sqrt(nansum(y .* y)) ;
	 ##w(I) = nan ;
	 ##w = Cy .* covest(w', w) ; 
	 w = Cy ;
      endif
      if any(isnan(w(:)))
	 error("too few y values") ;
      endif 

      try
	 Gy = chol(w) ;
	 I = true(1, columns(Gy)) ;
      catch
	 w = (w' + w)/2 ;
	 I = any(!isnan(w)) ; II = logical(I' * I) ;
	 w = reshape(w(II), sum(I), sum(I)) ;
	 Gy = chol(bend(w)) ;
      end_try_catch

      R = opp(x, y(:,I), GxI, Gy) ;
      wE(:,I) = nanmult(GxI, R, Gy) ;

   elseif (mode == "rgr")

      wE = nanmult(gpinv(x), y) ;

   endif 

   wE = nanmult(wE, Wgt^-1) ;

   if nargout > 1
      e = re(x*wE, y) ;
   endif 

   E(Jx,:) = wE ;

endfunction


function Q = opp (x, y, Gx, Gy)

   ## usage:  Q = opp (x, y, Gx, Gy)
   ##
   ## solve orthogonal procrustes problem

   [U S V] = svd(nanmult(Gy, y', x, Gx), "econ") ;
   Q = V * U' ;

endfunction


function y = gpinv (x, k)

   [m n] = size(x) ;

   if (nargin < 2), k = min(m, n) ; endif ;

   if (m > n)
      C = nanmult(x', x) ;
      [U S V] = svd(C) ;
      I = find(diag(S) > 0) ;
      I = I(1:length(I) <= k) ;
      U = U(:,I) ; S = S(I,I) ; V = V(I,:) ;
      CI = U * S^-1 * V' ;
      y = nanmult(CI, x') ;
   else
      C = nanmult(x, x') ;
      [U S V] = svd(C) ;
      I = find(diag(S) > 0) ;
      I = I(1:length(I) <= k) ;
      U = U(:,I) ; S = S(I,I) ; V = V(I,:) ;
      CI = U * S^-1 * V' ;
      y = nanmult(x', CI) ;
   endif 

endfunction


function c = covest (x)

   ## usage: c = covest (x)
   ##
   ##

   I = ~isnan(x) ;

   c = nanmult(x', x) ./ (I' * I) ;

endfunction
